real*8 function psi1(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi1=((r2-r)/(r2-r1))*((z2-z)/(z2-z1))
end
!*************************************************************************
real*8 function psi2(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi2=((r-r1)/(r2-r1))*((z2-z)/(z2-z1))
end
!*************************************************************************
real*8 function psi3(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi3=((r2-r)/(r2-r1))*((z-z1)/(z2-z1))
end
!*************************************************************************
real*8 function psi4(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi4=((r-r1)/(r2-r1))*((z-z1)/(z2-z1))
end
!*************************************************************************
real*8 function functionX(x1,x2,hx,x,i)
    real*8 x1,x2,hx,x
    k=mu_3d(i)
    if(k.eq.1) then
        functionX=(x2-x)/hx
    else
        functionX=(x-x1)/hx
    endif
end
!*************************************************************************
real*8 function functionY(y1,y2,hy,y,i)
    real*8 y1,y2,hy,y
    k=nu_3d(i)
    if(k.eq.1) then
        functionY=(y2-y)/hy
    else
        functionY=(y-y1)/hy
    endif
end
!*************************************************************************
real*8 function functionZ(z1,z2,hz,z,i)
    real*8 z1,z2,hz,z
    k=ku_3d(i)
    if(k.eq.1) then
        functionZ=(z2-z)/hz
    else
        functionZ=(z-z1)/hz
    endif
end
!*************************************************************************
subroutine calculateValueForNotTermNodes(v,igT,jgT,ggT,n,m,nTerm,nNotTerm)
    dimension v(n),ggT(m),igT(nTerm+1),jgT(m)
    real*8 v,ggT
        do i=1,nTerm
            v(nNotTerm+i) = 0.d0
            do j=igT(i),igT(i+1)-1
                v(nNotTerm+i)=v(nNotTerm+i)+v(jgT(j))*ggT(j)
            enddo
        enddo
    end
!*************************************************************************
subroutine mult_sq(g,b,v,n)
    dimension g(n),b(n,n),v(n)
    real*8 g,b,v

    do i=1,n
        g(i)=0.d0
        do j=1,n
            g(i)=g(i)+b(i,j)*v(j)
        enddo
    enddo
end    
!*************************************************************************
integer function searchElement2D(vertex,nVer,r,z,nodesCount,elemsCount)
    integer elemsCount
    real*8 vertex
    real*8 r,z,r1,r2,z1,z2
    dimension vertex(2,nodesCount),nVer(6,elemsCount)

    if(r.lt.1d-3) r=1d-3

    searchElement2D=1
    r1=vertex(1,nVer(1,searchElement2D))
    r2=vertex(1,nVer(2,searchElement2D))
    z1=vertex(2,nVer(1,searchElement2D))
    z2=vertex(2,nVer(3,searchElement2D))
    do while(searchElement2D+1.le.elemsCount.and.(r1.gt.r.or.r2.lt.r.or.z1.gt.z.or.z2.lt.z))
        searchElement2D=searchElement2D+1
        r1=vertex(1,nVer(1,searchElement2D))
        r2=vertex(1,nVer(2,searchElement2D))
        z1=vertex(2,nVer(1,searchElement2D))
        z2=vertex(2,nVer(3,searchElement2D))
    enddo
end
!*************************************************************************
integer function searchElement3D(vertex,edges,nodesforedges,x,y,z,nodesCount,&
    elemsCount,edgesCount)
    integer elemsCount,edgesCount,edges
    real*8 vertex
    real*8 x,y,z,x1,x2,y1,y2,z1,z2
    dimension vertex(3,nodesCount),edges(12,elemsCount),nodesforedges(3,edgesCount)

    searchElement3D=1
    jx=edges(1,searchElement3D)
    jy=edges(5,searchElement3D)
    jz=edges(9,searchElement3D)
    x1=vertex(1,nodesforedges(1,jx))
    x2=vertex(1,nodesforedges(2,jx))
    y1=vertex(2,nodesforedges(1,jy))
    y2=vertex(2,nodesforedges(2,jy))
    z1=vertex(3,nodesforedges(1,jz))
    z2=vertex(3,nodesforedges(2,jz))
    do while(searchElement3D+1.le.elemsCount.and.(x1.gt.x.or.x2.lt.x.or.&
                                                  y1.gt.y.or.y2.lt.y.or.&
                                                  z1.gt.z.or.z2.lt.z))
        searchElement3D=searchElement3D+1
        jx=edges(1,searchElement3D)
        jy=edges(5,searchElement3D)
        jz=edges(9,searchElement3D)
        x1=vertex(1,nodesforedges(1,jx))
        x2=vertex(1,nodesforedges(2,jx))
        y1=vertex(2,nodesforedges(1,jy))
        y2=vertex(2,nodesforedges(2,jy))
        z1=vertex(3,nodesforedges(1,jz))
        z2=vertex(3,nodesforedges(2,jz))
    enddo
end
!*************************************************************************
subroutine e2d(r,z,vertex2D,nVer2D,vRZ,vXYZ,dimV,nodesCount2D,elemsCount2D,nEl)
    integer elemsCount2D,dimV
    real*8 vertex2D,vertex3D,vRZ,vXYZ
    real*8 r,z,r1,r2,z1,z2,psi1,psi2,psi3,psi4
    dimension vertex2D(2,nodesCount2D),nVer2D(6,elemsCount2D),vRZ(nodesCount2D,dimV),vXYZ(dimV)

    r1=vertex2D(1,nVer2D(1,nEl))
    r2=vertex2D(1,nVer2D(2,nEl))
    z1=vertex2D(2,nVer2D(1,nEl))
    z2=vertex2D(2,nVer2D(3,nEl))

    do j=1,dimV
        vXYZ(j)=psi1(r1,r2,z1,z2,r,z)*vRZ(nVer2D(1,nEl),j)+&
                psi2(r1,r2,z1,z2,r,z)*vRZ(nVer2D(2,nEl),j)+&
                psi3(r1,r2,z1,z2,r,z)*vRZ(nVer2D(3,nEl),j)+&
                psi4(r1,r2,z1,z2,r,z)*vRZ(nVer2D(4,nEl),j)
    enddo
end
!*************************************************************************
subroutine inputParams_task3d(params,sig3d,nMaterials)
    integer sig3d
    real*8 params
    dimension params(2,*)

    nMaterials=0
    do while(.not.EOF(sig3d))
        read(sig3d,*)k,params(1,k),params(2,k)
        if(nMaterials<k) nMaterials = k
    enddo
end
!*************************************************************************
subroutine inputParams_task2d(params,sigma_txt,nMaterials2d)
    integer sigma_txt
    real*8 params
    dimension params(*)

    nMaterials2d = 0
    do while(.not.EOF(sigma_txt))
        read(sigma_txt,*)k,params(k)
        if(nMaterials2d<k) nMaterials2d = k
    enddo
end
!*************************************************************************
!*************************************************************************
!*************************************************************************